PTEN R Notebook

Multiplex assessment of protein variant abundance by massively parallel sequencing VAMP-seq - multiplex assay that uses fluorescent reporters to measure the steady-state abundance of protein variants in cultured human cells (each cellexpresses a single variant directly fused to EGFP…the stability of the variant dictates the abundance of the EGFP fusion and, accordingly, the green fluorescence signal of the cell) - used to assess PTEN and TPMT variants

# graph VAMP-seq scores relative to variant position in protein
#pten
pten1_hbond <- pten1_proc[!is.na(pten1_proc$hbond_sum),]
pten1_hbond$secondary_struct <- ifelse(is.na(pten1_hbond$helix), "unknown",
                        ifelse(pten1_hbond$helix==1, "helix",
                        ifelse(pten1_hbond$sheet==1, "sheet",
                        ifelse(pten1_hbond$helix==0, "neither",
                        "unknown"))))
pten_plot_hbond <- ggplot(pten1_hbond, aes(x=hbond_sum, y=score, colour=secondary_struct))+ geom_point(alpha=0.4) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("PTEN scores in relation to hydrogen bonding") + scale_color_manual(values=c("#FF4848", "#696969", "#5757FF")) + labs(colour="Secondary Structure") + theme_bw()
plot(pten_plot_hbond)

pten_plot_hbond1 <- ggplot(pten1_hbond, aes(x=hbond_sum, y=score))+ geom_point(alpha = 0.2) + ylab("VAMP-seq score")+xlab("DSSP Sum of hydrogen bonds")+ggtitle("PTEN scores in relation to hydrogen bonding") + theme_bw()
# was in aes, ggplot function call ---> colour=secondary_struct
#scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) + labs(colour="Secondary Structure")+
plot(pten_plot_hbond1)

#less hydrogen bonds ~ higher abundance

#method1 (basic, for visualizing in rstudio)
grid.newpage()
grid.draw(rbind(ggplotGrob(tpmt_dssp_schematic), ggplotGrob(tpmt_pos_mean), ggplotGrob(tpmt_heat), size = "last"))

Removed 1 rows containing missing values (geom_segment).Removed 1 rows containing missing values (geom_point).

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_dssp_schematic), ggplotGrob(pten_pos_mean), ggplotGrob(pten_heat), size = "last"))

Removed 12 rows containing missing values (geom_errorbar).

#plotting mean score vs variant changed to 
tpmt_end_sum <- summarySE(tpmt1_proc_wt, measurevar="score", groupvars="end")
tpmt_end_mean <- ggplot(tpmt_end_sum, aes(x=end, y=score)) +
  geom_bar(position=position_dodge(), stat="identity", colour="#999999") +
  geom_errorbar(aes(ymin=score-sd, ymax=score+sd), width=0.001, position=position_dodge()) +
  ylab("mean abundance") + xlab("variant amino acid") + theme_bw()
plot(tpmt_end_mean)

pten_a_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "A"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Alanine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "A"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()
pten_a_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "A"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Alanine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "A"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()
pten_r_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "R"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Arganine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "R"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()
pten_r_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "R"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Arganine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "R"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()
pten_n_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Asparagine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()
pten_n_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Asparagine variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()
pten_d_spread1 <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Aspartic Acid variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()
pten_d_aa <- ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Aspartic Acid variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()
pten_n_hydrodiff <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "N"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Asparagine variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "N"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()
pten_d_hydrodiff <- ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "D"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Aspartic Acid variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "D"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()
plot(pten_a_spread1)

plot(pten_a_aa)

plot(pten_r_spread1)

plot(pten_r_aa)

plot(pten_n_spread1)

plot(pten_n_hydrodiff)

plot(pten_n_aa)

plot(pten_d_spread1)

plot(pten_d_hydrodiff)

plot(pten_d_aa)

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Proline variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("Proline variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "P"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("Proline variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "P"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("T variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=end, colour = end)) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=end), scale = "width") + xlab("Amino acid") + ggtitle("T variant abundance scores")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=end, y=score), alpha = 0.75, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = (hydro2-hydro1))) + geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, start== "T"), aes(group=factor(position)), scale = "width") + xlab("Position in PTEN") + ggtitle("T variant abundance scores w/ hydrophobicity change")+ geom_point(data=subset(pten1_proc_wt, start== "T"), aes(x=factor(position), y=score), alpha = 0.85, position=jitter2) + scale_color_distiller(palette = "Spectral") + theme_bw()

#in order to find position with low stddev (excluding nonsense)
pten_variance <- summarySE((subset(pten1_data, end != "X")), measurevar="score", groupvars="position", na.rm=TRUE)

NaNs produced

#adjust N (minimum # variants at a position) and stddev
#5, 0.11
#10, 0.15
pten_variance_filtered <- subset(subset(pten_variance, N > 8 ), sd < 0.11)
pten_variance_filtered$position

 [1]  16  18  20  21  75 204 252 274 276 277 294 329 387 388 397

#nonsense variant scores are graphed as dots, but are not included in the overlaying violin plots
ggplot(no_nonsense, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered$position), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered$position), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(no_nonsense, aes(y=score, x=position, colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered$position), aes(group=position%%450), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered$position), aes(x=position, y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

#diff parameter 
pten_variance_filtered1 <- subset(subset(pten_variance, N > 10), sd > 0.25)
no_nonsense1 <- subset(pten1_proc_wt, end != "X")
#nonsense variant scores are graphed as dots, but are not included in the overlaying violin plots
ggplot(no_nonsense1, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense1, no_nonsense1$position %in% pten_variance_filtered1$position), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered1$position), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

ggplot(no_nonsense1, aes(y=score, x=position, colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense1, no_nonsense1$position %in% pten_variance_filtered1$position), aes(group=position%%450), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered1$position), aes(x=position, y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

chosen <- c(61, 68, 105, 108, 123, 127, 130, 132, 135, 155, 165, 173, 174, 246, 323, 335)
ggplot(pten1_proc_wt, aes(y=score, x=factor(position), colour = end)) + 
  geom_violin(draw_quantiles=c(0.5), data=subset(pten1_proc_wt, pten1_proc_wt$position %in% chosen), aes(group=factor(position)), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Clinvar path/likely path variant positions") +
  geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% chosen), aes(x=factor(position), y=score), alpha = 0.85, position=jitter1) + scale_color_manual(values=twenty_color) + theme_bw()

#pten_pos_mean <- ggplot(pten_sum, aes(x=position, y=score))+ geom_bar(position=position_dodge(), stat="identity", colour="#999999") + geom_errorbar(aes(ymin=score-sd, ymax = score+sd), width=1, size=0.3, position=position_dodge()) +ylab("VAMP-seq score")+theme(axis.title.x = element_blank(), axis.text.x = element_blank()) + scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) + geom_vline(xintercept=185, color="black", size=.1) + geom_vline(xintercept=350, color="black", size=.1)
#pretty trashy
ggplot(no_nonsense, aes(y=score, x=position)) + 
  geom_bar(data=subset(no_nonsense1, no_nonsense$position %in% pten_variance_filtered1$position), position=position_dodge(), stat="identity", colour="#999999") +
  geom_errorbar(data=subset(no_nonsense1, no_nonsense$position %in% pten_variance_filtered1$position), aes(ymin=score-sd, ymax=score+sd), width=1, size=0.3, position=position_dodge()) +
  #geom_violin(draw_quantiles=c(0.5), data=subset(no_nonsense1, no_nonsense1$position %in% pten_variance_filtered1$position), aes(group=position%%450), scale = "width") + 
  xlab("Position in PTEN") + ggtitle("Intolerant and tolerant amino acid positions") +
  #geom_point(data=subset(pten1_proc_wt, pten1_proc_wt$position %in% pten_variance_filtered1$position), aes(x=position, y=score), alpha = 0.85, position=jitter1) +
  #scale_color_manual(values=twenty_color) + 
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0))+theme(axis.title.x = element_blank(), axis.text.x = element_blank(), legend.position="none")

pten_variance_filtered2 <- subset(subset(pten_variance, N > 10), sd > 0.35)
ggplot(no_nonsense, aes(y=sd, x=position)) +
  geom_bar(data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered$position), position=position_dodge(), stat="identity", colour="#999999")

#all sd over .25, N>10
ggplot(pten_sum, aes(y=sd, x=position)) +
  geom_point(data=subset(pten_sum, pten_sum$position %in% pten_variance_filtered1$position)) +
  #geom_errorbar(data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered2$position), aes(ymin=score-sd, ymax = score+sd), width=1, size=0.3, position=position_dodge()) + 
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0))+theme(axis.title.x = element_blank(), axis.text.x = element_blank(), legend.position="none")

#this is how _ is defined: pten_variance <- summarySE((subset(pten1_data, end != "X")), measurevar="score", groupvars="position", na.rm=TRUE)
pten_variance_filtered4 <-subset(pten_variance, N>5)
variance_bar <- ggplot(pten_variance_filtered4, aes(y=sd, x=position)) +
  geom_segment(aes(x=position, xend=position, y=0, yend=sd), color="grey68") +
  geom_point(size=0.5) +
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0))+theme(axis.title.x = element_blank(), axis.text.x = element_blank(), legend.position="none")
plot(variance_bar)

#all sd
ggplot(pten_sum, aes(y=sd, x=position)) +
  geom_segment(aes(x=position, xend=position, y=0, yend=sd), color="grey48") +
  #geom_point() +
  #geom_errorbar(data=subset(no_nonsense, no_nonsense$position %in% pten_variance_filtered2$position), aes(ymin=score-sd, ymax = score+sd), width=1, size=0.3, position=position_dodge()) + 
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0))+theme(axis.title.x = element_blank(), axis.text.x = element_blank(), legend.position="none")

NA

pten_sum_filt <- pten_sum
pten_sum_filt$sd2 = pten_variance_filtered4[match(pten_sum_filt$position, pten_variance_filtered4$position), "sd"]
pten_pos_colored_mean <- ggplot(pten_sum_filt, aes(x=position, y=score))+ geom_bar(position=position_dodge(), stat="identity", colour=pten_sum_filt$sd2) + geom_errorbar(aes(ymin=score-sd, ymax = score+sd), width=1, size=0.3, position=position_dodge()) +ylab("VAMP-seq score")+theme(axis.title.x = element_blank(), axis.text.x = element_blank()) + scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) + geom_vline(xintercept=185, color="black", size=.1) + geom_vline(xintercept=350, color="black", size=.1)+ scale_colour_gradient(low = "white", high = "yellow")
plot(pten_pos_colored_mean)

#method2 (use for final layout, size specification, download)
gd=ggplot_gtable(ggplot_build(variance_bar))
maxWidth = grid::unit.pmax(ga$widths, gb$widths, gc$widths, gd$widths)
ga$widths <- as.list(maxWidth)
gb$widths <- as.list(maxWidth)
gc$widths <- as.list(maxWidth)
gd$widths <- as.list(maxWidth)
grid.newpage()

#storing, with specified widths!!
pdf('pten_tpmt_mean_heat_variance.pdf', width=8, height=6)
#grid.arrange(arrangeGrob(gC,gA,gB,nrow=3,heights=c(.1,.3,.8)))
grid.arrange(arrangeGrob(gc,gd,ga,gb,nrow=4,heights=c(.1,.15,.15,.6)))

dev.off()

quartz_off_screen 
                3 

#Identifying items in tail to investigate
pten1_nonsense <- subset(pten1_proc, class == "nonsense")
tpmt1_nonsense <- subset(tpmt1_proc, class == "nonsense")
pten1_synon <- subset(pten1_proc, class == "synonymous")
tpmt1_synon <- subset(tpmt1_proc, class == "synonymous")
pten1_no_missense <- subset(pten1_proc, class == "synonymous" | class == "nonsense")
ggplot(pten1_nonsense, aes(x=score)) + geom_histogram(binwidth=.01, colour="blue", fill="white") + theme_bw()

#+ geom_density()
ggplot(pten1_synon, aes(x=score)) + geom_histogram(binwidth=.01, colour="red", fill="white") + theme_bw()

ggplot(pten1_proc_wt, aes(x=score)) + geom_histogram(data=subset(pten1_proc_wt,class == "nonsense"), fill = "red", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_proc_wt,class == "synonymous"), fill = "blue", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_proc_wt,class == "missense"), fill = "green", alpha = 0.2, binwidth=.01) + theme_bw()

ggplot(pten1_no_missense, aes(x=score)) + geom_histogram(data=subset(pten1_no_missense,class == "nonsense"), fill = "red", alpha = 0.5, binwidth=.01) + geom_histogram(data=subset(pten1_no_missense,class == "synonymous"), fill = "blue", alpha = 0.5, binwidth=.01) + theme_bw()

ggplot(tpmt1_synon, aes(x=score)) + geom_histogram(binwidth=.01, colour="red", fill="white") + theme_bw()

ggplot(tpmt1_nonsense, aes(x=score)) + geom_histogram(binwidth=.01, colour="blue", fill="white") + theme_bw()

#0.55
nonsense_tail <- subset(pten1_nonsense, score > 0.6)
synon_tail <- subset(pten1_synon, score < 0.6)
nonsense_tail$secondary_struct <- ifelse(is.na(nonsense_tail$helix), "unknown",
                        ifelse(nonsense_tail$helix==1, "helix",
                        ifelse(nonsense_tail$sheet==1, "sheet",
                        ifelse(nonsense_tail$helix==0, "neither",
                        "unknown"))))
synon_tail$secondary_struct <- ifelse(is.na(synon_tail$helix), "unknown",
                        ifelse(synon_tail$helix==1, "helix",
                        ifelse(synon_tail$sheet==1, "sheet",
                        ifelse(synon_tail$helix==0, "neither",
                        "unknown"))))
#data[row,column]
n_tail <- nonsense_tail[,c(1,2,7,30,127)]
s_tail <- synon_tail[,c(1,2,7,30,127)]
n_tail$bp_pos <- (n_tail$position-1)*3
s_tail$bp_pos <- (s_tail$position-1)*3
n_tail

s_tail

#just in case there is a discernible pattern
s_tail_pos <- ggplot(s_tail, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Secondary Structure")+ggtitle("PTEN synonymous variant tail scores in relation to protein structure") + geom_vline(xintercept=47, color="black", size=.1) + geom_vline(xintercept=78, color="black", size=.1) + geom_vline(xintercept=122.5, color="black", size=.1) + geom_vline(xintercept=140, color="black", size=.1) + geom_vline(xintercept=165, color="black", size=.1) + geom_vline(xintercept=194, color="black", size=.1) + geom_vline(xintercept=209, color="black", size=.1) + theme_bw()
plot(s_tail_pos)

#help visualizing NMD rules
n_tail_pos <- ggplot(n_tail, aes(x=position, y=score, colour=secondary_struct))+ geom_point(size=.3) + scale_x_continuous(minor_breaks = seq(0, 405, 5)) + scale_color_manual(values=c("#FF4848", "#00C853", "#5757FF", "#A9A9A9")) +ylab("VAMP-seq score")+xlab("Position in PTEN")+labs(colour="Secondary Structure")+ggtitle("PTEN nonsense variant tail scores in relation to protein structure") + geom_vline(xintercept=47, color="black", size=.1) + geom_vline(xintercept=78, color="black", size=.1) + geom_vline(xintercept=122.5, color="black", size=.1) + geom_vline(xintercept=140, color="black", size=.1) + geom_vline(xintercept=165, color="black", size=.1) + geom_vline(xintercept=194, color="black", size=.1) + geom_vline(xintercept=209, color="black", size=.1) + theme_bw()
plot(n_tail_pos)

TPMT_abun_CADD <- ggplot(tpmt_merge, aes(x=abundance_class, y=CADD_raw_rankscore)) + geom_violin(draw_quantiles = c( 0.5))+ylab("CADD raw rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_CADD)

TPMT_abun_SIFT_conv <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(SIFT_converted_rankscore))) + geom_violin(draw_quantiles = c(0.5))+ylab("SIFT conv rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_SIFT_conv)

TPMT_abun_POLY <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(Polyphen2_HDIV_rankscore))) + geom_violin(draw_quantiles = c( 0.5))+ylab("Polyphen2 HDIV rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_POLY)

TPMT_abun_POLY1 <- ggplot(tpmt_merge, aes(x=abundance_class, y=as.numeric(Polyphen2_HVAR_rankscore))) + geom_violin(draw_quantiles = c( 0.5))+ylab("Polyphen2 HVAR rankscore")+xlab("Abundance Class") + theme_bw()
plot(TPMT_abun_POLY1)

Pred_abun_SIFT <- ggplot(tpmt_merge, aes(abundance_class)) + geom_bar(aes(fill = SIFT_pred)) + ggtitle("Abundance class vs SIFT prediction of Damaging or Tolerated") + theme_bw()
plot(Pred_abun_SIFT)

trial_sep <- tpmt_merge[c(21,23,24,26)]
tpmt_merge_expand <- separate_rows(tpmt_merge, c("Polyphen2_HDIV_score", "Polyphen2_HDIV_pred", "Polyphen2_HVAR_score", "Polyphen2_HVAR_pred"))
Pred_abun_HVAR <- ggplot(tpmt_merge_expand, aes(abundance_class)) + geom_bar(aes(fill = Polyphen2_HVAR_pred)) + ggtitle("Abundance class vs Polyphen2 HVAR predictions") + labs(subtitle = "D: Probably Damaging, P: Possibly Damaging, B: Benign") + theme_bw()
plot(Pred_abun_HVAR)

twenty_color1 = c("#D02028", "#A4C33B","#53958B", "#E6A3B4", "#C5A0CA", "#554DA0", "#99247E", "#402059", "#82421B", "#7E807E", 'black', "#EDD941", "#F2F08E", "#EEC898", "#E1A12F", "#76C158",  "#BCDDAE", "#85782E", "#315935", "#A1DAE0", "#486EB6")
pten_dssp_schematic1 <- ggplot() +
  geom_segment(aes(x = 1, y = 0, xend = max(pten_extra$position)), yend = 0, size = 1, color = "grey70") +
  geom_point(data = subset(pten_extra, !is.na(xca)), aes(x = position, y = 0), color = "black", size = 1.8) +
  geom_point(data = subset(pten_extra, sheet == 1), aes(x = position, y = 0), color = "pink", size = 1.5) +
  geom_point(data = subset(pten_extra, helix == 1), aes(x = position, y = 0), color = "cyan", size = 1.5) +
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) +
  scale_y_continuous(breaks = NULL, expand = c(0,0)) +xlab("Position in PTEN") + ylab("\n \n \n") +
  theme(panel.border = element_blank(), axis.text.y = element_blank())

aas <- c("A", "C", "P", "X")
aas1 <- c("S", "C", "P", "X")
aas2 <- c("A", "C", "D", "E", "F", "G", "H", "I", "K", "L", "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
pten_a_pos <- ggplot(data=subset(pten1_proc_wt, end=="A"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="A" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Alanine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
#plot(pten_a_pos)

pten_s_pos <- ggplot(data=subset(pten1_proc_wt, end=="S"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="S" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Serine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)

pten_c_pos <- ggplot(data=subset(pten1_proc_wt, end=="C"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="C" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Cysteine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)


pten_x_pos <- ggplot(data=subset(pten1_proc_wt, end=="X"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="X" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of nonsense variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_a_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_s_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_c_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_x_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

#snake at the end (arches over grey region)
pten_p_pos <- ggplot(data=subset(pten1_proc_wt, end=="P"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="P" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Proline variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_p_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

twenty_color1 = c("#D02028", "#A4C33B","#53958B", "#E6A3B4", "#C5A0CA", "#554DA0", "#99247E", "#402059", "#82421B", "#7E807E", 'black', "#EDD941", "#F2F08E", "#EEC898", "#E1A12F", "#76C158",  "#BCDDAE", "#85782E", "#315935", "#A1DAE0", "#486EB6")
pten_dssp_schematic1 <- ggplot() +
  geom_segment(aes(x = 1, y = 0, xend = max(pten_extra$position)), yend = 0, size = 1, color = "grey70") +
  geom_point(data = subset(pten_extra, !is.na(xca)), aes(x = position, y = 0), color = "black", size = 1.8) +
  geom_point(data = subset(pten_extra, sheet == 1), aes(x = position, y = 0), color = "pink", size = 1.5) +
  geom_point(data = subset(pten_extra, helix == 1), aes(x = position, y = 0), color = "cyan", size = 1.5) +
  scale_x_continuous(breaks = seq(0, 403, 20), expand = c(0,0)) +
  scale_y_continuous(breaks = NULL, expand = c(0,0)) +xlab("Position in PTEN") + ylab("\n \n \n") +
  theme(panel.border = element_blank(), axis.text.y = element_blank())
aas <- c("A", "C", "P", "X")
aas1 <- c("S", "C", "P", "X")
aas2 <- c("A", "C", "D", "E", "F", "G", "H", "I", "K", "L", "M", "N", "P", "Q", "R", "S", "T", "V", "W", "Y")
pten_a_pos <- ggplot(data=subset(pten1_proc_wt, end=="A"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="A" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Alanine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
#plot(pten_a_pos)
pten_s_pos <- ggplot(data=subset(pten1_proc_wt, end=="S"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="S" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Serine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
pten_c_pos <- ggplot(data=subset(pten1_proc_wt, end=="C"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="C" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Cysteine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
pten_x_pos <- ggplot(data=subset(pten1_proc_wt, end=="X"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="X" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of nonsense variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_a_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_s_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_c_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

grid.newpage()
grid.draw(rbind(ggplotGrob(pten_x_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

#snake at the end (arches over grey region)
pten_p_pos <- ggplot(data=subset(pten1_proc_wt, end=="P"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="P" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Proline variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_p_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

#matching snake at the end!
pten_g_pos <- ggplot(data=subset(pten1_proc_wt, end=="G"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="G" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Glycine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_g_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

pten_h_pos <- ggplot(data=subset(pten1_proc_wt, end=="H"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="H" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Histidine variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_h_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))

#2nd from last grey region
pten_w_pos <- ggplot(data=subset(pten1_proc_wt, end=="W"), aes(x=position, y=score, colour=start))+ geom_point(data=subset(pten1_proc_wt, end=="W" & start %in% aas2), size=.6) + scale_x_continuous(breaks=seq(0, 403, 20), minor_breaks = seq(0, 403, 5), limits = c(0, 403), expand = c(0,0)) + scale_y_continuous(expand = c(0,0), limits = c(-0.25, 1.5)) + scale_color_manual(values=twenty_color1) +ylab("Abundance score")+labs(colour="Reference amino acid")+xlab(NULL)+ggtitle("Abundance scores of Typtophan variants") + theme(axis.text.x = element_blank(), legend.position='top') + geom_hline(yintercept=1, color="black", size=.1)
grid.newpage()
grid.draw(rbind(ggplotGrob(pten_w_pos), ggplotGrob(pten_dssp_schematic1), size = "last"))